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SUMMARY 


A robust, iterative solution p>rocedure is presented for the parabolized 
Navier-Stokes or higher-order bound ary- 1 aver equations as applied to <ub- 
sonic viscous-inviscid interaction flows. The robustness of the present pro- 
cedure is due. in part, to an improved algorithmic formulation. The present 
formulation is based on a re-interpretation of stability requirements for 
this class o'* algorithms and requires only second-order-accurate backward, 
or central differences for ail streamwise derivatives. I'pstream influence 
is provided for through the algorithmic formulation and iterative sweeps 
in j. The primary contribution to robustness, however, is the boundary- 
condition treatment, which imposes global constraints to control the con- 
vergence path. Discussed are successful calculations of subsonic, strong 
viscous-inviscid interactions, including separation. These results are con- 
sistent with Navier-Stokes solutions and triple-deck theory. 


SYMBOLS 


C/ skin friction coefficient, 2Tjpeul 
Cp pressure coefficient 

e total energ}’ per unit mass, m^/sec^ 

rfi mass-flow rate, kg/sec 

p static pressure, kg/m-sec^ 

R radius of curvature, m 

Rol Reynolds number, uLju 

u velocity, ^-component, m/sec 

V velocity, ^/-component, m/sec 

X Cartesian coordinate, m 

y Cartesian coordinate, m 

6* displacement thickness, m 

f) transformed coordinate, f? = ?y(j, y) 

p djTiamic viscosity, kg/m sec 

1/ kinematic viscosity, m^/sec 

^ transformed coordinate, ^ = x 

p density, kg/m"^ 

r shear stress, pdujdy 

'I' streamfuDction, / pu 

Subscripts: 

CL wake centerline value 

e edge value 

i grid index in x 

j grid index in y 

TE trailing-edge ^'alue 

■w wall value 

Superscript: 

n iteration level index 
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INTRODUCTION 


AerodjTiamic flows often involve physical phenomena not properly ac- 
counted for in first-order boundary-layer theory (ref. 1), even when coupled 
to an outer inviscid solution. For strong viscous-inviscid regions^ such as 
the trailing-edge of an airfoil, or in a shock/boundary-layer interaction, 
streamline curvature leads to a significant normal pressure gradient dp/dy 
within the boundary layer. Additionally, upstream efiTects may become im- 
portant even within the thin viscous shear layer. Separation may also occur, 
further complicating the solution within these flow regions. 

Some coupled methods (refs. 2-4) modify the results of their boundary- 
layer solutions within such regions in an attempt to correct for these deficien- 
cies of first-order boundary-layer theory. However, the appropiateness of 
these corrections, \*alid only to first order, clouds related issues, such as tur- 
bulence modeling, within these regions. To resort to the time-dependent 
Navier-Stokes approach for these regions seems a needlessly expensive way 
of accounting for the weak ellipticity of these flows. Problems with grid 
resolution are also encountered in the application of time-dependent Navier- 
Stokes methods owing to storage and CFL condition limitations, even for 
the implicit methods prevalent (ref. 5). The present work addresses this 
clear need for an eflBcient viscous-solution method, one that is intermediate 
between currently available first-order bound ary-layer and Navier-Stokes 
techniques, is capable of treating upstream influence and small-scale separa- 
tion, and is suitable for viscous-inviscid coupling. 

Previous workers have formulated sophisticated schemes with various 
regimes of applicability to solve viscous layers with higher-order effects. 
Here we concentrate on spatially marching schemes because of their implied 
promise of minimizing both storage requirements and convergence times 
while retaining the capability of treating complex geometries. Parabolized 
Navier-Stokes (PNS) methods have demonstrated considerable success w hen 
applied to supersonic flow over complex bodies. The most prevalent PNS 
methods (refs. 6,7) currently are single-sweep, noniterative and are faster 
than time-dependent Navier-Stokes solvers by orders of magnitude, thus 
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making three-dimensioDal computations quite practical. Such PNS methods 
effectively solve the continuity, i-momentum, y-momentum, and energ}^ 
equations for p, pu, pv, and pe, respectively. Pressure is obtained from the 
equation of state. Currently, such PNS solvers are restricted to supersonic 
flows with thin subsonic regions. The treatment of the pressure terms within 
the subsonic regions is only approximate but is critical to obtain stable 
marching. For stable marching, a minimum streamwise step Ax > Ax^in 
must be taken. This minimum step size is roughly the height of the subsonic 
region. Axial separation is not allowed. These PNS methods currently in- 
volve only single sweeps in the x-direction and thus cannot provide for 
upstream influence. 

Partially-parabolized Navier-Stokes (PPNS) methods preside solutions 
for large regions of subsonic and even reversed flow. In these methods (refs. 
8,9), the X- and y- momentum equations are spatially marched with the 
pressure field assumed known. The continuity equation is then used as a 
basis for corrections to the pressure field. Some PPNS methods also solve 
an elliptic Poisson equation which aids in transmitting pressure information 
throughout the entire flow field, but also adds to the computational time 
and effort. The specification of u and ? boundary conditions does not easily 
fit into current coupling schemes. Upstream influence is possible with these 
methods and multiple sweeps in the x-direction are necessary. 

More suited to application in a viscous-inviscid coupling scheme are 
PNS methods that are essentially boundary-layer methods with incorpora- 
tion of the t/-momentum equation. Incorporation of the y-momentum equa- 
tion elevates the conventional first-order boundary-layer equations into the 
PNS set of equations. Upstream influence is provided through multiple- 
sweeps in the x-direction and through proper design of the algorithm. 
Primary contributions to these methods were made by Rubin (ref. 10), 
Lin and Rubin (ref. 11), Maghoub and Bradshaw (ref. 12), and Chen and 
Bradshaw (ref. 13). These PNS methods effectively solve the continuity, 
x-momentum, y-momentum and energv' equations for v, u, p, and pe respec- 
tively. The equation of state is used to solve for p. For incompressible flows, 
the variable pe and the energj' equation may be dropped. 

Rubin applied such a PNS/boundary-layer method to an incompres- 
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sible separation problem. The method is implicit with the equations coupled. 
Rubin emphasized forward-differencing of dp/dx inconsistent with the differ- 
encing treatment of the remaining streamwise derivatives as a requirement 
for stable marching. Extension to compressible flows was accomplished by 
Lin and Rubin in a study of supersonic viscous flow over an inclined cone. 
This class of flows contains negligible upstream influence and is efficiently 
solved by the single- sweep PNS methods. Convergence of their multiple- 
sweep PNS/boundary-layer method did occur, albeit slowly. Rubin and 
Reddy (ref. 14) subsequently applied a multigrid technique to aid conver- 
gence and robustness on fine grids. 

Chen and Bradshaw presented an efficient viscous-inviscid coupling 
scheme based on the PNS/boundary-layer method of Maghoub and Brad- 
shaw which includes the t/-momentum equation. Good agreement of the 
calculations with experimental Cp versus x was obtained for two airfoils 
at transonic speeds but at a small angle of attack. The boundary-layer 
equations are solved uncoupled, and it would appear that the technique, 
being based on the method of characteristics, cannot calculate through 
flow-reversal. Chen and Bradshaw concede that if d‘^pfdxdy is significant 
their method may fail. 

Both of the above higher-order boundary-layer methods are significant 
contributions. Both of these methods, however, have limitations as to ac- 
curacy and robustness. Accuracy is limited in both methods since forward- 
differencing of dp/dx is used in the j-momentum equation for stability 
and upstream influence, and backward-differencing is used for all other x- 
derivaiives. As a model equation, even the differential form of the Bernoulli 
equation will be solved inaccuratel}'. 

In this paper, the requirement of forward-differencing of dpfdx for 
stability and upstream influence is removed. A consistent algorithm is for- 
mulated. thereby improving the inherent accuracy of this class of methods. 
Coarser grids are allowed which aids the convergence rate, and the robust- 
ness of the current method is considerably enhanched relative to the prior 
methods, a result of both the dp/dx and boundary-condition treatments. 

The method in this paper is an implicit, finite-difference solution of Ih 
PN.S or higher-order boundary-layer approximation of the Navier-Stokes 


equations. In addition to the continuity and j-momentum equations solved 
by first-order boundary-layer techniques, inclusion of the ^/-momentum equa- 
tion is accomplished. Inclusion of the y-momentum equation allows for solu- 
tion with normal pressure gradients but converts the parabolic first-order 
boundary-layer equations into an elliptic set. This is demonstrated by drop- 
ping the viscous terms, thereby obtaining the Euler equations as a subset. 
.\n otherwise parabolic spatially marching procedure is rendered elliptic by 
providing for upstream propagation of information relating to the pressure 
field. Second-order- accurate backward or central differences are used for 
the streamwise derivatives, leading to an accurate and consistent scheme. 
Physical insight allows the x- and j/- momentum equations to be written in a 
form so that the scheme is also stable. Iterative sweeps in the 2 -direction are 
required with global storage for only pressure being neccessary. The only 
special treatment required for separation is to set pudujdx = 0 in the x~ 
momentum equation. 

The boundary conditions of the method presented here are similar 
to those of conventional first-order boundary- layer methods with either a 
direct edge-pressure, or an inverse mode. EflScient algorithms exist for 
coupling with an inviscid code, using either the direct or the inverse boun- 
dary conditions. The boundary conditions implemented are thus well-suited 
for incorporation of the method into a coupled viscou.s-inviscid procedure. 
The psuedodirect boundary condition is actually an adaption of a proposed 
semi-inverse coupling procedure. For duct flow, where mass flow remains 
constant, the present code can proceed to a solution without the necces- 
sity of coupling. The influence of the type of boundary condition used 
on convergence is observed and discussed. To validate the method, results 
are presented for a variety of flow conditions. Comparisons are made with 
available Navier-Stokes solutions and triple-deck theory. 
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NUMERICAL FORMULATION 


The equations that are solved for incompressible, t'wo-dimensional flows 


are 


du dv 
dx dy 


= 0 


du , du 

cpu— -|-pv— = 

dx dy 

dv , dv 

dx dy 


dp 

dx 

?L 

dy 



( 1 ) 


where the FL.ARE approximation (ref. 15) is used 


c=l if u > 0 
= 0 if u < 0 


Boundary conditions are discussed later. 

For coupled systems of equations, a t'ariable that is being solved for 
may be associated with a particular equation through diagonal dominance of 
the resultant matrix system. As with conventional first-order boundary- 
layer techniques, the solution for u and i is associated with the j-momentum 
and continuity equations, respectively. The solution for pressure i<5 as- 
sociated with the ^-momentum equation. In contrast, PPNS methods solve 
for u and ? from the x- and y- momentum equations, respectively. This 
difl'erence is not just a subtle point. Not only are simpler boundary condi 
tion specifications required in the present formulation, but also the pressure 
field is being solved in an implicit coupled fashion simultaneously with the 
velocity field. In PPNS methods, the pressure-field calculations are de- 
coupled and are lagged relative to the velocity-field calculations. Generally 
speaking, methods that are implicit and fully coupled are more robust and 
converge faster than methods that are not. 

One key idea in the present formulation is that in deriving the finite- 
difference representation of the above differential equations, no equation is 
written ai the same location. By so doing, stable marching is achieved while 
also arriving at a sparse block-tridiagoiial matrix system w hich may be more 
economically solved. Figure 1 depicts the various nodes involved in writing 
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the finite-difference equations. In figure 1, the equation is written formall}" 
at {i,j-\/ 2 ) for continuity; (i,j) for ^--momentum; and (t-l,y + 1/2) for 
j/-m omentum. In marching the solution one streamwise step, the variables 
u and r are being solved for at 2,, and p is being solved for at 2,-1. The 
continuity and 2-momentum equations use backward differences for the 2- 
derivatives. The backward-difference expression for the dp/dx term with 
evenly-spaced .42 is 


ap _ P7-2,j 

lAx 


( 2 ) 


The superscript n refers to the global iteration level. Note that p”-2,j 
has been solved for, p” j ^ is being simultaneously solved for from the y- 
moinentum equation, and p^'J must be obtained from a previous global 
iteration. Upstream values of at 2, are involved in the y-momentum 
solution for p” at 2,-1. Upstream influence is thus provided for through the 
(t,<-momentum and dpjdx) pressure formulation. W^re the y-momentum 
equation to be vritten at 2,, not only would the upstream propagation of 
pressure information be suppressed, but the scheme would be unstable for 
small Ax. 

The metric treatment used allows more general body shapes than 
does a rectangular grid, while minimizing the complexity involved during 
the early phases of code development. Consider the rather simple metric 
transformation 

x = y = yU,ri) 

^ = 2 ; p==Ti{x,y) = j 
The 2- and y- derivatives may be written as 

d _ dr] d d 

dx dx dr] d^ , 

dy dy dr] 


The term dr]fdx can be written for the above transformation as 


■which leads to 


d d dyd dy t/t,j — yt-\,j 

The 2 -derivatives in the differential system, equation 1, are then written 
using equation 6, and the y-derivatives are written as y-derivatives using 
higher-order-accurate numerical forms for variable grid spacing. The solu- 
tion then actually occurs in a hybrid y system. Metric treatments of the 
tj-pe presented by Steger (ref. 16) are contemplated for future incorpora- 
tion; however, a quite wide range of problems with interesting geometries 
may be studied with the metric treatment described. 

using the above metrics, second-order- accurate finite differences, and 
the finite-difference molecule depicted in figure 1, the differential syst.em 
given by equation 1 can be written in matrix form 


)y( )j* 1 




(T) 


where 

j —ij'max refers to the jth differencing moIecule(see fig.l) 
k refers to the equation; 

= 1, continuity 
=2, 2 -momentum 
=3, y-momentum 
m refers to the ^'ariable: 

= 1 ) ~ ~ 

= 2 , {U2)j = 

=3, (?73)j = Pfl.3 

The interior elements for matrices A,B,C, and D are given in the appendix. 
This matrix system is then solved using a standard 3x3 block-tridiagonal 
matrix solver rewritten to take advantage of the matrix sparsity. The 
sparsity of the matrix system is depicted in figure 2. The convective 
coefficients, and {pv)x,j, appearing in the i-momentum equation 

involve unknowns being solved for. Iteration for these coefficients makes 
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the finite-difference approximation to the differential system second-order- 
accurate. Typically, three iterations are used with the values at (2-1,^) 
used in the first iteration. Extrapolation for these coefBcients was tried; 
however, the converged solutions for separated flows exhibited weak spatial 
oscillations. Note that the j/-momentum equation is solved for at {i-l.j + 
1 /2); therefore, these convective coefficients are known. 

A computer code written in FORTRAN implements the above al- 
gorithm for two-dimensional steady, incompressible laminar flows. The 
code operates in either a PNS or a first-order boundary-layer (BL, where 
dp/dy = 0 replaces the ^/-momentum equation) option. Such an ability 
to easily choose between the PNS and BL formulations allows accurate 
et^luation of the significance of the ^-momentum equation for various flow 
conditions. Several boundary-condition options are implemented. These 
are next described. 

The boundary conditions for a solid surface are 

continuity at Xj : = 0 at wall 

j-momentum at x. : u, i = 0 at wall ,,, , 

1 ' ( 8 | 
on 

— = 0 at edge 
dy 

An additional edge condition is required for the y-momentum equation 
on the pressure. The treatment of this additional edge condition leads 
to four boundary-condition options: the direct mode, duct mode, inverse 
mode, and psuedodirect mode. 

Direct mode - The direct mode applies an input pressure distribution 
Pt,jmax versus X, as an edge condition to the y-momentum equation. This 
option is referred to as the PNS-Direct or the BL-Direct mode depending 
on whether the PNS or the BL option is in effect. 

Duct mode - The duct mode applies a correction to the overall pressure 
level at x, to maintain a constant mass flow, m, through the duct. The value 
of m is obtained from the inlet profile at X]. With P'l-\,jmax' 

assumed known, the following FORTR.AN assignment statements 
are used: 
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( 9 ) 



( used as the j/-momentum equation edge condition) 


The pressure level correction, PC”, is found by iteration through 
a modified secant method relating to PC”. Typically three or four 
iterations are required. This option is referred to as the PNS-Duct or BL- 
Duct mode. 

Inverse mode - The inverse mode is similar to the duct mode except 
that the pressure correction being applied is adjusted so that an input edge 
streamfunction is obtained, versus where 


This option is referred to as the PNS-CS’A" or BL-E\A’’ mode. 

Psuedodirect mode - The psuedodirect mode is similar to the inverse 
mode except that the edge streamfunction distribution versus 2 :, is 
adjusted every global iteration. The adjustment is accomplished to give, 
upon convergence, a desired edge- velocity distribution; it is based on a 
coupling algorithm that is being explored. The algorithm used is 



This inverse mode is similar to the more conventional (^’-specified inverse 
mode. Note that 'I'g and 6 * are related by 


'Pe = «e(ye “ Vxv ~ 


an 


q,n ^ ^n-1 _f_ 


-AUe 


( 12 ) 


where 



( 13 ) 
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Since 66* fbue is not easily evaluated, this term is at present dropped and 
nnderrelaxation is used giving 




n-l 


^ ^n-1 , 

c e » 


U 


(H) 


Here Au^ is the difference between the desired edge velocity at a:, 
and the edge velocity from the last global iteration. The desired 
edge velocity can be either an input edge- velocity distribution, as with the 
present psuedodirect mode, or can be an edge- velocity distribution obtained 
from an inviscid solver, as in a viscous-inviscid coupling procedure. A 
similarity exists to Carter’s 6* semi-inverse algorithm (ref. 17) 


(5 ” = 6 




c, inviscid 


(15) 


The psuedodirect option demonstrates the essential features of the 
coupling algorithm and accelerates convergence over the direct mode for 
the specified edge-velocity boundary condition. The psuedodirect mode is 
used only with the PNS option and is referred to as the PNS-PSD mode. 
The logic for the PNS-PSD mode is demonstrated in block form in figure 
3. 


For the symmetric wake, the boundary condition for the 2 -momentum 
equation differs at y = 0. The symmetry of u about y = 0 is made use of 
to evaluate the j-momentum equation at y = 0 as 


/ ,n I n in 19, 


(>'"2- ""1) 


(16) 


The d/dx terms are expanded as before. The symmetric wake boundary 
condition above is second-order accurate and gives a smoother represen- 
tation of the inner wake than does the first-order du/dy = 0 condition 
normally applied. 
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RESULTS 


The method described in the previous sections was applied to several 
laminar flows that involve strong viscous-inviscid interactions. These flows 
included the Briley laminar separation cases and the symmetric laminar 
wake of a flat plate. Much of the development work on the algorithm 
and boundary-condition treatments was accomplished while considering the 
Briley separation cases. 


Laminar Separation Results 


Briley (ref. 18) solved the time-dependent Navier-Stokes equations over 
a flat plate for a set of four linearly retarded flows. The edge velocity is 
prescribed at yc = 3.81 x 10“® m for these four flows as 


Ue(m/sec) = 30.48(1 — 3x), x < Xc 
= C , X > Xc 


(171 


The four Briley cases difiTer in the location Xc and velocity level C. Brile) 
faired the curve at Xc to avoid a discontinuity in slope. The kinematic vis- 
cosity is L' = 1.4864 X lO""^ m^/sec. The specification of these parameters 
by Briley was in English units (the nondimensional Reynolds number being 
the important constant) and have been subsequently converted to metric 
for this study. Two of the four Briley cases (Briley-1 and -2) are attached 
flows; the other two (Briley-3 and -4) exhibit separation bubbles of differing 
extent. Elliptic or upstream effects are significant in these flows and, in 
contrast, a pressure-prescribed first-order boundary-layer solution (Howarth 
flow ) predicts separation for all four cases. Since the geometry of these cases 
is simple, and separation with upstream influence does occur, these flows 
have been used in various studies as standard test cases for validating and 
developing new numerical algorithms for calculations in fluid mechanics. 

The initial efforts at solving the Briley cases made use of an algorithm 
quite similar to that of Rubin, in that forward-differencing of dpfdx and 
backward-differencing of dujdx were implemented in the i-momentum 
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equations. Coarse-grid solutions inevitably did not separate since the ad- 
verse pressure gradient was being relieved prematurely owing to the use of 
a forward-differenced dpfdx. Fine-grid solutions could be obtained with 
underrelaxation but the method lacked robustness. It was realized during 
these initial efforts that for these methods to be generally and routinely 
applicable for industrial and aerodjuamic design purposes, developments 
had to be made to 

1. Eliminate the inconsistent use of forward-differenced dpjdx and 
baclrward-differenced dujdx in the i-momentum equation 

2. Enhance accuracy on coarse grids 

3. Enhance robustness, particularly when applied to flows with severe 
adverse pressure gradients 

4. Promote rapid convergence 

An important development was the realization, based partly on physi- 
cal arguments, that stable marching of this class of algorithms requires that 
the pressure 'V’ariable solved for be the trailing pressure (with a dominant 
negative coefficient) in the dpjdx term; Rubin’s stability requirement of 
forward-differencing for the dpjdx term is unneccessarily restrictive. Solving 
for the trailing pressure leads to damping of error terms, whereas solving 
for the leading pressure leads to error amplification and the associated 
departure solutions described by Light hill (ref. 19). Sub5:equently, the 
present fully implicit algorithm was formulated with second-order- accurate 
backward-differenced d jdx used exclusively for the continuity and the x- 
momentum equations, and second-order-accurate central-differenced d jdx 
used for the y-momentum equation. Improvements in accuracy and robust- 
ness were demonstrated by the elimination of the forward-differenced dpjdx. 
The convergence rate was, however, still felt to be slower for the Briley flows 
than one might reasonably expect. The development of a plot package 
depicting streamline and pressure contours revealed intermediate solutions 
with substantial streamline cun^ature, which in some cases verged on being 
physically unrealistic. At this point, a decision was made to develop a treat- 
ment of the boundary conditions that would more affirmatively control the 
convergence path of these intermediate solutions. This led to the PNS-ESrV’' 
and later to the PNS-PSD boundary-condition modes. 
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The Briley-3 separated- flow case was selected for more intensive study 
instead of the Briley-4 case because of limitations of the PNS-Direct mode. 
AJso, 'Various Navier-Stokes methods give somewhat different solutions for 
the Briley-4 case since significant viscous effects extend to the edge of 
the computational domain, and details of boundary-condition treatment 
acquire importance. Fully converged solutions for the Briley-4 case were 
obtained with the PNS-PSD mode, however. Thus, the attached Briley- 1 
and the separated Briley-3 cases were chosen for use in the present stud}' 
for validation purposes. Five modes of operation of the present code were 
applied to the Briley-3 and -1 cases: (1) PNS-Direct; (2) PNS-ES"^'; (3) 
PNS-PSD; (4) BL-Direct; and (5) BL-IN^^ Coarse- and fine-grid solutions 
were obtained. The Briley-3 case is discussed first. 

Briley-S case - Figure 4 depicts the edge- velocity distribution used for 
the Briley-3 case as the input boundary condition for the PNS-Direct and 
PNS-PSD modes. The position for separation as given by the Briley Navier- 
Stokes and by the first-order boundar}'-layer (Howarth fiow) solutions are 
marked; reattachment as calculated by Briley is indicated as well. Murphy 
(ref. 20) demonstrated that the difference between the Navier-Stokes and 
first-order boundary-layer separation points is a function of the ReA'Uolds 
number of the flow, with agreement occurring at sufficiently high Reynolds 
number. For the PNS-tNV and BL-EN^’’ modes the input boundary con- 
dition used was the edge streamfunction calculated from the PNS-PSL' 
mode. The output result for the inverse modes is a calculated edge- velocity 
distribution. The resultant edge- velocity distribution from the PNS-ENA’ 
mode is identical with the input distribution to the PNS-PSD mode, as 
one might reasonably expect. The BL-INA' mode, however, gives an edge- 
velocity distribution in which the greatest differences occur in the region of 
the separation bubble, where dpfdy ^ 0 effects are most significant. 

The streamline and pressure contours obtained for the PNS-INA' mode 
are combined in figure 5. The PNS-PSD solution is essentially identical. 
The streamlines are not evenly spaced, being more tightly spaced close to 
the surface to provide more detail. The pressure contours are evenly spaced 
in pressure. The j-spacing of the pressure contours gives an indication of 
the streamwise pressure gradient, and the slope of a particular pressure 
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contour is inversely related to the normal pressure gradient: 


_ ^P, / dy 


(18) 


A fine-grid (36 i-stations) solution is shown. Minor difiTerences in stream- 
lines for the coarse-grid (18 2 -stations) solution exist, but are not significant. 
For most purposes, the accuracy of the coarse-grid 10th iteration solu- 
tion would suffice. Thus, rapid convergence to an accurate solution is 
demonstrated on a coarse grid for the current algorithm with the PNS-E\~V^ 
boundary-condition mode. 

The combined streamline and pressure-contour plot makes it possible 
to gain insight into both the flow being computed and the algorithm employ- 
ed. One connection between the pressure field and streamlines is through 
the y-momentum equation, rewritten here as 

= _ djvlu) 

dy R ’ R dx ^ 

where R is interpreted as the radius of curvature of the local streamline. 
Were R to be interpreted as the surface curvature, as in some boundary- 
layer methods, the ellipticity of the present method would be suppressed. 
Note that close to the surface within 6*, even though the streamline cur- 
vature is substantial, the normal pressure gradient is nearly zero. This is, 
of course, because the pu^ term is quite small. Only in the inviscid region 
and in the viscous region outside <5* do strong dpfdy effects occur, being 
most significant for those 2 -stations in the vicinity of the separation bubble, 
where streamline curvature is large. In the vicinity of the separation point, 
the streamline curvature and associated dp/dy effects tend to relieve the 
surface-pressure gradient {dpjdx) relative to the edge-pressure gradient, 
thereby delaying separation relative to a first-order boundary-layer solution. 
The weak adverse pressure gradient is, however, prolonged at the surface, 
thereby delating reattachment till well past the 2 -station where the edge- 
pressure gradient is relieved. It is interesting to note that reattachmeut 
occurs in the presence of a weak adverse pressure gradient at the surface 
Reattachment is apparently accomplished by the shear-stress transport of 
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momentum from the outer regions predominating over the weak residual 
adverse pressure gradient. 

The considerable advantage in convergence rate of the PNS-EVA' mode 
over that of the PNS-Direct mode can best be understood in terms of this 
interaction between the pressure field and the streamline contours. With the 
edge streamline specified, the streamlines within the computational domain 
approximate fairly well the streamlines for the final converged solution, even 
on the first global iteration. The dp/dy field establishes itself rather quickly 
in the global iteration process. For the PNS-Direct mode, the streamlines 
within the computational domain are not so constrained, and, particularly 
for separated flows, the initial streamlines differ markedly from the final 
converged solution; there is slower establishment of the pressure field, with 
slower overall convergence. 

Figure 6 shows velocity profiles u versus y obtained using the PNS- 
PSD mode. In the interest of clarity, not all ^-locations are presented. It is 
obvious that the velocity levels in the reversed-flow region are quite small, 
thereby justifjing the FLARE approximation used in the present study. It 
is also apparent that the viscous region approaches but does not quite reach 
the edge of the computational domain for this Briley-3 case. ,\Jso observe 
that the outflow profile is Blasius in shape, in agreement with Briley. As 
indicated in figure b, a stretched y-spacing is used, with 14 of 40 nodes 
adequate to resolve the inlet boundary layer. Solutions with even spacing 
in y give nearly the same result. 

Figure 7 presents the distribution of skin friction Cf versus i and dis- 
placement thickness 6* versus x calculated by the PNS-EsW (identical to the 
PNS-PSD solution) and the BL-LNW modes. Also given are the correspond- 
ing results from Briley’s calculations. The skin-friction results indicate that 
the PNS-Cs"V mode and Briley calculations agree identically with respect to 
reattachment location, with the PNS-LNY mode giving separation one-half 
grid point farther downstream. The BL-ESW mode gives surprisingly good 
agreement with Briley’s separation point, with reattachment occcurring two 
grid points upstream of the Briley calculation results. Good agreement for 
the 8* distribution exists for the PNS-ESW and the Briley solutions. The 
BL-ESW method once again demonstrates some differences for Cf and (5* 
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in the region of significant dpjdy. The BL-INV method used appears to 
be a particularly attractive calculation method in viev? of the single global 
sweep required. Fletcher (ref. 21) has reached similar conclusions concern- 
ing an inverse BL, ^’-prescribed, finite- difference method. An error in edge 
velocity in regions of significant dp/dy inevitably occurs, however, with 
these inverse BL methods. 

For the Briley-3 case, the PNS-PSD and PNS-IN\^ modes give results 
quite comparable to the Navier-Stokes solutions by Briley. These methods 
converge rapidly, with about 20 global iterations required for the fine-grid 
PNS-ENT solution and 10 global iterations for the coarse-grid PNS-IN\’ 
solution. The iterations required for the PNS-PSD mode depend on the 
initial distribution. Using an initial distribution deduced from the 
Briley results, the fine-grid PNS-PSD converges in 40 global iterations. 
Here convergence is indicated when no significant changes occur in C/. 
5* , or Uf. The pressure residual for these two methods goes to essentialh 
machine zero. The PN.S-Direcl mode was also applied to the Briley-3 case 
with results similar to those of the PNS-PSD mode. The convergence 
characteristics are not as desirable, with the maximum pressure residual 
leveling off before reaching machine zero. For separated flows, the PNS- 
PSD mode exhibits considerable advantage over the PNS-Direct mode with 
respect to the convergence rate and robustness. 

Briley-1 case - The Brilej-l attached-flow case was also calculated bt 
the PNS and BL methods. Figure 8 indicates the edge- velocity distribution 
Uf versus i specified by Briley for this case. This Ue distribution is similar to 
that for the Briley-3 case but the adverse pressure gradient is not sustained 
for as long. Navier-Stokes solutions do not separate, even though separa- 
tion is indicated by first-order boundary-layer theory (Howarth solution). 
Thus, upstream influence is significant for this flow. The skin-friction and 
displacement thickness results are shown in figure 9. The PNS-Direct, PNS- 
LV\', and PNS-PSD modes all give results in agreement with the Navier- 
Stokes solutions of Briley. Once again, the BL-INN^ mode gives quite accep- 
table results with the greatest discrepancy occurring in regions of significant 
dpjdy. .Also indicated are BL-Direct solutions for this case in the regioE 
leading to separation. The separation point indicated by the BL-Direci 
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results is in good agreement with the Howarth separation point. It is in- 
teresting to note that with the FLARE approximation being used, the BL- 
Direct mode will march through separation and reattachment. The solu- 
tions past separation are highly nonphysical, however, with a tremendous 
ejection of material upward from the separation point (see the discussion of 
Landau and Lifshitz (ref. 22) regarding the singularity at separation of the 
first-order boundary-layer equations). The substantial streamline curvature 
in the BL-Direct solutions close to separation is at considerable variance 
with the first-order bound ary- layer theory assumption of dp/dy = 0. 


Trailing Edge Results 


The laminar, flat-plate trailing-edge provides a further example of 
the applicability of the present PNS formulation to strong viscous-inviscid 
interactions. In the absence of upstream influence, a simple discontinuity 
in boundary condition occurs for first-order boundary-layer theory at the 
trailing-edge. Goldstein (ref. 23) provides such first-order boundary-layer 
solutions for the trailing-edge/near-wake problem. However, for the Navier- 
Stokes equations and equation subsets which properly account for upstream 
influence, the laminar trailing-edge problem contains a singularity at the 
trailing-edge. Stewartson (refs. 24,25) points out, based on the triple-deck 
subset of the Navier-Stokes equations, that the skin friction has an inverse 
square-root singularity as the trailing-edge is approached. The displacement 
thickness will, however, vary smoothly through the trailing-edge region. 
Interacting boundary-layer theory, where upstream influence is accounted 
for by the boundary-condition formulation, appears to adequately calculate 
this flow. Melnik et al. (ref. 2) has demonstrated a similar theory for a 
turbulent trailing edge. To better understand the trailing-edge problem 
and as a test case involving strong viscous-inviscid interaction, calculations 
using the PNS and BL formulations described earlier are presented below. 

The computation is for a fiat plate, 1 m long, with Re^ = 10""^^. The 
computational domain is for the interval between x = 0.5 and x = 1.5 
m, with the trailing edge located at i = 1. Two values for the edge of 
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the computational domain are considered. Both of these remain outside 
the "viscous layer. The inlet conditions are Blasius, but with 8* adjusted 
to account for the leading-edge interaction. The grid is unevenly spaced 
with 59 2 -stations and 40 j/-nodes, and is concentrated in the t railing- 
edge region where Ax = 0.01. Veldman (ref. 26) has calculated the 
above problem earlier with both interacting boundary-layer and triple-deck 
solution techniques. 

Several approaches are applied to this problem. These include 

1. PNS-PSD, with j/e = 0.12 m and Ug = 100 m/sec 

2. BL-Direct, with surface pressures from PNS-PSD(^e = 0.12) 

3. BL-Direct, with dp/dx = 0 and Ug = 100 m/sec; this case rep- 
resents noninteracting, first-order boundary-layer theory with a combined 
Blasius flat-plate and Goldstein near-wake solution. 

Initial PNS calculations were made with the computational edge at 
= 0.03, just outside the viscous region. .Although the solution did exhibit 
upstream influence, this tendency appeared damped with unacceptable 
quality of agreement with the triple-deck solutions. Two possibilities were 
explored in efforts to explain this disagreement. One involved the incorpora- 
tion of the remaining viscous terms into the algorithm(e.g., (fi'Uj)x, {pVy)y. 
and Stable spatial marching occurred for the full Navier-Stokes 

equations and the solutions actually presented are full Navier-Stokes solu- 
tions. However, these solutions showed only very minor changes from the 
PNS solutions and did not explain the disagreement above. With this ex- 
tension to the full Navier-Stokes equations, the organization of the present 
primitive- variable formulation is similar to the streamfunction/vorticity for- 
mulation of Murphy. 

The reason behind this disagreement was found, however, to lie in 
the second possibility. Examination of the pressure-contour plot for the 
PNS-PSD {pe = 0.03j case (see fig. 10) revealed pressure contours that 
appeared to be constrained by the edge of the computational domain. To 
alleviate this constraint, the computational edge was moved outward in 
several stages to yg = 0.12. The extent of upstream influence increased 
t.oward the triple-deck solution, after which no further changes occurred. 
Figure 11 depicts the pressure contours for the PNS-PSD (t/^ = 0.12) case 
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and should be contrasted with the data in figure 10. Note the different 
y- sc ales in figures 10 and 11. We see that the dist ortion of the pressure 
field by the trailing-edge propogates nearly to the inlet boundary for the 
ye = 0.12 case. 

Except for a small range of ^-stations very close to the trailing edge, 
much of the pressure-field adjustment occurs in the inviscid region and 
dpfdy 0 within the viscous region. The relative success of interacting 
boundary-layer theory for the trailing-edge problem is traceable to the 
validity of the assumption of a zero normal pressure gradient for much of 
the boundar}'-la3^er and wake viscous regions. To verify this, the surface 
pressures from the PNS-PSD (yg = 0.12) calculations were applied as the 
input to the BL-Direct mode. As will be seen, reasonable agreement with 
Veldman’s results are obtained. 

Figures 12 and 13 show distributions of the skin-friction coeflBcient 
Cj versus i, and the displacement thickness b* versus x for the three 
cases. Figure 14 shows the distribution of the normalized wake centerline 
velocity ucil'^t versus x. Also included is the triple-deck-based solution of 
Veldman. The PNS-PSD (y^ = 0.12) case shows excellent agreement with 
Veldman’s results, confirming that the present formulation is capable of 
correctly accounting for both the viscous and inviscid regions of this strong 
viscous-inviscid interaction flow. The results for the BL-Direct (with PNS- 
PSD surface pressure) case prove to be quite acceptable, with differences 
occurring for b* closest to the trailing edge, where the assumption of a 
zero normal pressure gradient is invalid. The skin-friction agreement is 
excellent, probably because the correct surface pressures are being applied. 
For the trailing-edge flow, the most significant effect of the viscous-inviscid 
interaction is to force the inviscid region to adjust the pressure field as 
applied to the viscous region. 

The source of this pressure-field readjustment is indicated in a com- 
parison of figures 15 and 16, which depict the streamline pattern for the 
noninteracting BL-Direct mode {dpjdx = 0) and the BL-Direct mode (with 
PN.S-PSD surface pressure). The streamline patterns for the PNS-PSD 
calculations are quite similar to those in figure 16. For the BL-Direct 
{dpjdx = 0) calculations, a discontinuity develops in the streamlines at 
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the trailing-edge. This streamline discontinuity is the result of the sudden 
acceleration of the innermost portion of the viscous layer upon release of the 
no-slip condition, a discontinuous change in boundary condition. Consider 
that the ^/-momentum equation can be written as dp/dy = —pu^/R; R 
is the radius of curvature of the local streamline. At the trailing edge, the 
discontinuity in streamlines implies that 0 and dp/dy oc. Thus, 
the streamline pattern generated by the BL-Direct {dp/dx = 0) calcula- 
tion is in sharp contradiction with the first-order bound ary- layer theory 
assumption that dp/dy = 0. If the j/-momentum equation is incorporated 
or the edge-pressure boundary condition is allowed to adjust to account for 
upstream infiuence, then the pressure field will adapt to smooth out the 
streamline curA'ature until a consistent solution is obtained. 


CONCLUDING REMARKS 


Progress on the development of a parabolized Navier-Stokes solution 
method suitable for application to strong viscous-inviscid interactions has 
been described. .An improved algorithm and boundary-condition treatment, 
which enhances rapid convergence to accurate solutions on coarse grids 
while retaining stable marching for small marching steps, is reported. To 
provide for proper upstream infiuence, the algorithmic formulation performs 
iterative sweeps in the j-direction and requires global storage only for pres- 
sure. The method has been ’V’ulidated through the successful computation 
of two subsonic, strong viscous-inviscid interactions, including a separation 
bubble and a trailing edge. Detailed comparison with Navier-Stokes solu- 
tions and triple-deck theory is made with excellent agreement indicated. 

The improved algorithm stems from a conceptual reinterpretation of 
stability requirements presented earlier by Rubin. In the present work, this 
stability requirement for forward-differencing of dp/dx is removed. Rather, 
this paper demonstrates that stable, subsonic marching of the PNS equa- 
tions merely requires that the pressure variable being solved for should be 
the trailing pressure in the dp/dx term. This idea is used to formulate the 
present algorithm; within each equation, consistent differencing is imple- 
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mented for all streamwise deri'V'atives. 

Boundary condition inodes for external flows are direct pressure, in- 
verse edge-streamfunction, and psuedodirect modes. With these boundary 
conditions, the present code is well suited for coupling with an inviscid 
code. For internal flows, global mass conservation is used. The boundary- 
condition treatment proves crucial in enhancing rapid convergence. For 
these strong viscous-inviscid interactions, the global constraints imposed 
by the inverse and psuedodirect boundary conditions provide the improve- 
ment in robustness over the direct boundary condition required to make 
this formulation a practical and useful approach. 

The computations for strong viscous-inviscid interactions demonstrates 
that the present PNS algorithm and boundary-condition formulation is 
a rapid, accurate, and robust approach which incorporates the essential 
features of elliptic upstream influence. For the trailing-edge flow, the PNS 
algorithm was extended to include the full Navier-Stokes terms with stable 
spatial marching achieved. To this author’s knowledge, this is the first 
successful spatially marched calculations of the full Navier-Stokes equations 
in primitive variable form for subsonic, strong viscous-inviscid interactions. 

The inclusion of a zero- and a two-equation turbulence model, and the 
incorporation of the compressible terms into the present code, recommen- 
dations for future work, have already been partially implemented. The 
combination of a turbulent compressible version of the present PNS code 
with an Euler solver in a coupling procedure should then be accomplished. 
The code would then serve as an efficient computational tool in turbulence 
modeling studies of strong viscous-inviscid interactions. 
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APPENDK: MATRIX ELEMENTS 


The matrix system being solved is of the form 

{■^k “1“ i^k,m)j{Um)j “f“ = (^k)j (-^1) 

Tfhere 

y ^Ifjmax refers to the jth differencing molecule 
k refers to the equation: 

=1, continuity 
=2, a;-momentum ' 

=3, J^momentum 
m refers to the variable: 

. — {^l)j — ~ 

= 2 , {U2)j = A'Uij = Uij , , 

=3, {Us)j = Pi-1, j 

The interior elements for the continuity equation are 

AVj-i : = 1 , 

AMj-1 : (Ai2)j = 0.5 • BXQ • DYX - DEXM 
Pi-ij-i : (Ai3)j = 0 


Avj 

(Bn)j 

= 1 


Auj 


= 0.5 • BXO 

DYX -f DEXM 

Pi-l,j 

(B.s), 

= 0 


Avj^i 

(Cn),; 

= 0 


AUy+1 

(C'l2)j 

= 0 


P»-i,y+i 

(Cl3); 

= 0 



(D,)j 





-K-u- 

-u^.,j.,)^DEXM 


+ 0.5 ■ DYX '- BX2 ■ 

- W2J + <-2,i-l)] 


(A2) 
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The interior elements for the a:-momentum equation are 


Avj-i 

Auj.i 

Pi-1, ri 



Pi-hj ' 

£,Vj+, : 
Auj+i : 
pi-1, j-\-l • 


(A2i)j = 0 

{A22)j = JRX • CYM - DTM 

{A23.)j — 0 
(B2i), = 0 

{B22)j = ip^)i,j • BXO + RX . CYO 
-\-{DTP-^DTM) 

(^s)y = BXl 

{C2l)j = 0 

{C22)j = RX • CYP - DTP 
(C2z)j = 0 


(r»2)i = -{P<‘)h ■ BX2 ■ 

+ KTP 
: + OTA/ 

- OTX ■ (cr p • pj'-Vi + cro ■ p?;/ 
+ cyA/-p?,-.l,) 
-(BAro.p”-' + BX 2 p?. 2 ,,) 


Note that the {pu)2j and RX terms are found through iteration. 



The interior elements for the y-momentum equation are 


: (^i)j = 0 

Auj-i : (A32)i = 0 

: (.433)y = 0 

Avy : = -0.5 . DYM • CXP • {pu)lij 

Atiy: (B32)/= 0.5^DYM‘CXP-ipv)l,j 

' (B33)y = — 1 

At;y+1 :• (C3i)> = ^0.5 DYM- CXP • {pu)f.^j^^ 

Auy+i : (C32)i = 0.5 • DYM • CXP • (ptO^.y+i 

Pt-i,y+i : (<^33)y= 1 

(Z) 2 )y = -0.5 . DYM • CXM • [(pi0r.i,y('^”2.y ~ <i.:) 

+ (pt')”i,y+i(t^r-2,y+i - ^r-i.y+i>] 
+ 0.5 . DYM . CXM • 

+ (p«)?-i,y+i{^r-2.y+i - 
+ 0.5 . DEXP . I(pt;)r.i,y + (P^;);i,y+i] 

• K-i.,^1 - <ij) 

- 0.5 . DEXP . [{pu)l,j + (pt;)” 



where 


DX = Xi— Xi-i 

DX2R — {xi — Xi.i)/{xi — Xi-2) 

BX\ = ~{Xi — Xi.2)l{Xi-i — Xi-2) 

BXO = -5X1/(1 - DX2R • DXIR) 

BX2 = -(BXO-f BXl) 

CXP = (Zi-i — Xi.2)l\{Xi — — x,-i)] 

CXM = — (a-,- — zm)/I( 2 ,- - Xi. 2 ){xi.i — a:, -. 2 )] 

CA"0 == -{CXP + CXM) 

P = , ivij y»,jf-l)/((2/j,; + l 2/»j'-0] 

CYM= ivij+i - yi,j)/[{yij - yij~i){yi,j+i - j/t,;-])] 
CYO = -{CYP-i-CYM) 

DYX = {yij yi,j-i)f{xi 

DYM = (y,-ij4-i - yi.ij) 

DEX = (y»,i 

DEXM = -0.5 • [(y,- j - 

2/t-i j'i)]/(2^i 

DEX P = 0.5 • {{yij j/tj+i) 

(y»-2,j 2/t-2,j'-i)]/(2’t 2 :^- 2 ) 

DTP = PA' • 

/[(y»ii+i yt,y-i) • (yi.y+i yt.y)] 

/{(y».y+i yi,j-i) ■ (y»,y y»,y-i^] 

KX - {pv)l - . PA + {pu)l^ . PPA^ 
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Figure 1. - Computational molecule. 
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Figure 2. - Matrix structure. 
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Figure 3. - Block logic for PNS-PSD mode. 
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Figure 4. - Edge- velocity distribution for Briley-3 case. 
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Figure 5. - Streamline and pressure contours for Briley-3 case: 
PNS-I^rV'^ solution. 
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F igure 6. - Velocity profiles for Briley-3 case: PNS-ENA" solution. 
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Figure 7. - Skin-friction and displacement-thickness distributions 
for Briley-3 case. 
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Figure 8. - Edge-velocity distribution for Briley-1 case. 
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Figure 9. - Skin-friction and displacement-thickness distributions 
for Briley- 1 case. 
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Figure 10. - Pressure contours for trailing-edge interaction: 
PNS-PSD with ye = 0.03 m. 
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Figure 11. - Pressure contours for trailing-edge interaction: 
PNS-PSD with ye = 0.12 m. 
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Figure 12. - Skin-friction coeflBcient distributions for trailing-edge 
interaction. 
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Figure 13. - DisplacemeDt-thickDess distributions for trailing-edge 
interaction. 
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Figure 14. - Wake centerline velocity distributions for trailing-edge 
interaction. 


44 



X, m 


Figure 15. - Streamline contours for trailing-edge interaction: 
BL-Direct “with dpjdj = 0. 
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Figure 16 . - Streamline contours for trailing-edge interaction: 
BL-Direct with PNS ■wall pressures. 


46 


-V- 


1. Report No. 2. Government Accession No. 

NASA TM-84378 

3. Recipient's Catalog No. 

4. Title and Subtitle 

PARABOLIZED NAVIER-STOKES SOLUTIONS OF SEPARATION 
AND TRAILING-EDGE FLOWS 

5. Report Date 

June 1983 

6. Performing Organization Code 

7. Author(s) 

James L. Brown 

8. Performing Organization Report No. 

A-9394 

10. Work Unit No. 

T-4219 

9. Performing Organization Name and Address 

NASA Ames Research Center 
Moffett Field, Calif. 94035 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

506-51-11-04-00 

16. Supplementary Notes 


Point of Contact: James L. Brown, Ames Research Center, MS 229-1, 

Moffett Field, Calif. 94035. (415) 965-6192 or FIS 448-6192. 


16 Abst'-act 

A robust, iterative solution procedure is presented for the parabolized 
Navier-S tokes or higher-order boundary-layer equations as applied to sub- 
sonic viscous-inviscid interaction flows. The robustness of the present 
procedure is due, in part, to an improved algorithmic formulation. The 
present formulation is based on a reinterpretation of stability requirements 
for this class of algorithms and requires only second-order-accurate back- 
ward or central differences for all streamwise derivatives. Upstream influ- 
ence is provided for through the algorithmic formulation and iterative 
sweeps in x. The primary contribution to robustness, however, is the 
boundary-condition treatment, which imposes global constraints to control 
the convergence path. Discussed are successful calculations of subsonic, 
strong viscous-inviscid interactions, including separation. These results 
are consistent with Navier-Stokes solutions and triple-deck theory. 
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